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Abstract 

Q ' We present the results of a numerical study of the 2 xL spin- 2 Heisenberg 

p\J . ladder. Ground state energies and the singlet-triplet energy gaps for 4 < L < 

14 and J±/ J\\ = 1 were obtained in a Lanczos calculation and checked against 

earlier calculations by Barnes et al. (even L < 12). A related moments 

^ ' technique is then employed to evaluate spin response functions for L = 12 

^ ' and a range of J±/J\\ (0 - 5). We comment on two issues, the need for 

reorthogonalization and the rate of convergence, that affect the numerical 

utility of the moments treatment of response functions. 



PACS number(s): 75.10.Jm,75.40.Gb,75.40.Mg 



Typeset using REVT^ 



Heisenberg spin ladders have attracted considerable attention recently [|l| due to possible 
connections to materials exhibiting high T^ superconductivity: theoretical studies 0] hint 
at the possibility of even-chain ladders becoming superconductors when doped with charge 
carriers. There is some experimental support for this possibility as Sro.4Ca13.6Cu24O41.84, a 
material with spin-| chains and 2-chain ladders, was shown to superconduct at 12 K and 3 
GPai. 

Spin ladders are also fascinating theoretically because of their unexpected behavior when 
viewed as interpolators between the spin-| ID antiferromagnetic Heisenberg chain and the 
2D square analog. The latter is fully ordered at low temperatures |^, while Bethe 
demonstrated in the 1930's that spin-spin correlations in the ID chain have a slow power- 
law decay. Yet the transition between these limits by forming spin-^ ladders with increasing 
numbers of legs is not smooth: ladders with even numbers of legs have a finite gap to the 
lowest triplet state and an exponential decay of spin-spin correlations, while odd-leg ladders 
have gapless excitations and a power-law fall off of spin-spin correlations 0. These issues 
and examples of materials exhibiting these properties are discussed in several recent reviews 
(see, e.g., ^). 

In the hope of gaining deeper insight into such systems, numerical modelers have em- 
ployed a variety of techniques to study spin ladders, including exact diagonalizations with 
the Lanczos algorithm [^ , quantum Monte Carlo simulations [^,0 , and approximate density 
matrix renormalization group methods using rung or plaquette bases (see White et al. @ 
and Piekarewicz and Shepard [§). The exact calculations, while limited to small L, play an 
important role in testing approximation schemes, and also in evaluating dynamic quantities, 
such as spin responses, that are difficult to treat in other approaches. In this report we 
present Lanczos results for the ground-state energy and singlet-triplet gap for 2 x L systems 
through L = 14, which we compare to the even L < 12 calculations [0] of Barnes et al. 
We then study the evolution of the dynamical spin response function for L = 12 as the 
ratio of rung to leg interaction strengths J_l/J\\ is varied. As this response can be mea- 
sured in inelastic neutron scattering experiments 0, it provides an important test of spin 
interactions in ladder materials. This response has also been evaluated in recent plaquette 
basis approximation schemes [§ . The present calculations are based on a Lanczos moments 
expansion that can be iterated to arbitrary accuracy. We discuss some numerical aspects of 
this procedure that are relevant to spin ladder calculations. 

The Hamiltonian for the spin-i Heisenberg spin ladder consisting of two coupled chains 



IS 



H=J±T. S^■S, + J\\ Y. S^■SJ, (1) 

where i is a lattice site on which one electron sits, {i,j)± denotes nearest neighbor sites on 
the same rung, and {i,j)\\ denotes nearest neighbors on either leg of the ladder. We used 
periodic boundary conditions along the legs of the ladder. The ratio J±/J\\, the relative 
strength of the rung and leg interactions, depends on the choice of material being modeled. 
In the strong coupling limit, J±/J\\ -^ 00, the electrons across each rung form an S* = 
pair, and the ground state wave function is the resulting product [0 . The ground state thus 
has S = and an energy/spin proportional to J±, with perturbative corrections of relative 
size J\\/J±- 



In the Lanczos algorithm |T^ the Hainiltonian is written in tridiagonal form recursively, 
using a series of operations 

H\vi) = /?i_i|fi_i) + ai\vi) + Pi\vi+i) (2) 

in which the next basis vector jt'i+i) is generated from the previous one, \vi), with the choice 
of |f i) depending on the application. We truncated this series after k steps and diagonalized 
the resulting k x k matrix by the QL algorithm |T^. The resulting energy per spin for the 



singlet ground state and gap to the first triplet state are compared to the results of Barnes 
et al. in Tables I and II. The agreement is very good, with only very minor differences 
appearing for large L. Our calculations through L = 13 were done in both single and double 
precision: the results are identical to the accuracy employed in the tables. The L = 14 
calculations were performed in single precision only. 

The ground-state energy per spin can be extrapolated to the bulk limit using a scaling 
function similar to that of Barnes et al. W\ 
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-L/Lo 



/(L)-/(oo) = Co(-l)^-^^ (3) 

where p = 2. A fit to the results of Table I yields Cq = 1.12, Lq = 3.86, and f{oo) = -0.5780. 
The even-L results for the spin-triplet gap extrapolate with p = 1 to f{oo) = 0.502 (Co = 
3.61, Lq = 3.82). 

We now consider the evolution of the dynamic spin response function S{q,uj), where q 
and uj are the three-momentum and energy transfer, as J±/J\\ is varied. S{q,uj) is defined 

by 

S{q, uj)=Y: I {n\S{q} | (?...) p5(a; - u^) (4) 

n 

where \g.s.) denotes the ground state and where the sum is taken over a complete set of 
excited states \n) of energy u;„. The spin transition operator is 

S{q) = Y.S,e^'-'', (5) 

where the sum extends over all sites. In particular, if the dynamic spin response is probed 
at g = (g^., qy) = (tt, tt) = q^^T,, in units of the inverse lattice spacing, then 

s{q.n) = j:i-^ysj, (6) 

so that the operator sign alternates from site to site. 

In the Barnes et al. work Lanczos techniques were employed in evaluating S{qT,.„,uj) for 
L=8 [0. Quite recently Piekarewicz and Shepard |Q studied the L = 6, 8, and 16 systems 
by exploiting a plaquette truncation of the basis. Thus one motivation for the present effort 
is to provide a series of exact calculations in somewhat larger systems (L=12) that can serve 
as benchmarks for approximate methods, like those of Piekarewicz and Shepard, that are 
now being applied to dynamical quantities. 



The Lanczos method is particularly well suited to the evaluation of inclusive response 
functions. Once an initial Lanczos expansion has been carried out to the point where the 
ground state is fully converged, the vector 

Sziq7r7T)\9-S.) (7) 

can be formed and its norm determined. The resulting normalized vector \vi) can then be 
used as the starting vector in a second Lanczos expansion, which is then stopped after k 
iterations. If one denotes the resulting eigenvectors and eigenvalues of the fc-dimensional 
Lanczos matrix by |/j) and e^, then 

k 

3Y.mSMn)\9-s.)\'S{u;-e,), (8) 

viewed as a distribution in ut, reproduces the lowest 2/c — 1 moments of the exact distribution 
given in Eq. (4) [0. Thus the broad outline of the response function is determined after a 
few iterations, with finer details emerging as the addition of higher moments increases the 
resolution. The Lanczos moments technique for inclusive response functions is thus exact in 
two senses: the lowest 2k — 1 moments are correctly determined, and for any specified limit 
of resolution (e.g., that achieved in some experiment) the iterations can be continued until 
a sufficient number of moments are obtained to produce an overall profile that is exact at 
the scale. 

Figs, la and lb give the dynamic spin response function per spin, S{qT,T,,u!)/2L for L = 
12, J±/J\\ = 1, and J\\ = 1, smoothed by a Gaussian resolution function with a standard 
deviation of 0.05. These initial calculations were done to determine, for this choice of 
resolution, the required number of iterations. These results, and those of Fig. 2, suggest 
that ~ 70 iterations are needed to produce a fully converged distribution. This conclusion, 
however, depends on one's choice of J_l/J\\'- when this ratio is increased, the proportion of 
the response carried by high-lying excitations drops, while at the same time the distribution 
of strength at high u exhibits more structure. Both effects slow the rate of convergence. 
Thus we found it necessary to use 200 iterations in the case of J±/J\\ = 5.0. 

A second numerical issue is the absence of exact orthogonality of the Lanczos vectors 
when Eq. (2) is implemented numerically [|rT|. Errors associated with the overlaps of a newly 
generated Lanczos vector \vi) with previous vectors can be quite troublesome: spurious over- 
laps with extremum eigenvectors can grow, in successive iterations, because they contribute 
so strongly to higher moments. This can lead, for example, to repeated reconvergence of the 
ground state and to distortions in inclusive response functions. In many applications this 
difficulty makes repeated reorthogonalization by the Gramm-Schmidt procedure necessary, 
a step that becomes costly when a large number of iterations are performed. However, the 
need for reorthogonalization varies greatly from application to application. The results from 
our exploration of this issue for Heisenberg spin ladders is shown in Fig. 2, where a. k = 
70 calculation with reorthogonalization in each iteration is compared to fc = 70 and 500 
calculations without. The calculations were performed for L = 10 and J±/J\\ = 1. The 
Heisenberg ladder Hamiltonian appears to be remarkably immune to numerical orthogonal- 
ity difficulties: no differences among the three calculations are readily discernable. Thus the 
remainder of the calculations reported here were done without a reorthogonalization step. 



Table III and Figs. 3a-l give our main results, S^q-^^T,, u) per spin for J\\ = 1.0 and J±/J\\ 
ranging from 0.0 to 5.0. The distributions have been smoothed by a Gaussian resolution 
function with a = 0.05. As the momentum transfer corresponds to the inverse lattice size, 
the operator reverses the orientation of nearest neighbor spins: a low-lying spin triplet state 
increasingly dominates the spin response function as J± is increased (see Table III). The 
gap between the singlet ground state and the strong triplet state increases with increasing 
J±/J\\, in agreement with the strong coupling (large J±) prediction of Egap ~ J± — J\\. 

The strength above the first triplet state is always modest, starting at ~ 16 % for J±/J\\ 
= and declining montonically to ~ 0.02 % for J±/J\\ = 5. The pattern of this strength, 
however, becomes more distinctive with increasing J±/J\\. Thus in principle this part of 
the dynamic spin response, while accounting for little of the total strength, could be used 
in combination with the singlet-triplet gap to test whether real materials respond as simple 
spin ladders. 

The total response strength per spin is not a monotonic function of J±/J\\, but instead 
increases from J±/J\\= to a peak at about J±/J\\ ~ 0.5, then declines steadily above 
this value. The results are shown in Table III. We also examined the evolution of the 
total strength, and strength carried by the first triplet state, as a function of L for fixed 
J±/J\\ = 1-0. The fraction of strength carried by the first triplet state appears to converge 
rapidly: the resuhs for L= 6, 8, 10, and 12 are 97.66%, 97.04%, 96.81%, and 96.73%, 
which extrapolate to a bulk limit of 96.70%. The bulk limit of the absolute strength per 
spin carried by the first triplet state is somewhat less certain, as this quantity appears to 
be rather slowly approaching an asymptote in our calculations: the results for L = 6, 8, 
10, and 12 are 3.26745, 3.52202, 3.66170, and 3.73619. These values can be fit very well 
by a scaling function of the form of Eq. (3) with p = 0, leading to a bulk limit of 3.83. 
To compare with |^, we repeated these calculations for open boundary conditions, finding 
2.80696, 2.99800, 3.12141, and 3.20260. The L = 8 resuh agrees with the corresponding 
exact calculation of |Q, who also used their approximate plaquette basis calculations to 
estimate a bulk limit open-boundary-condition value of 3.30-3.36. (Note that the values 
given in ||^ have been multiplied by 3 to take into account the different normalization of 
S{qT^T,,uj) used there.) If we employ the same scaling function in this case, our calculations 
yield a bulk limit 3.35, compatible with this range. Unfortunately we have results for too 
few L values to test numerically whether our assumed scaling function is reasonable. Thus 
there could be substantial errors in these estimates. 

Finally, we stress that the dynamical spin response calculations described here are quite 
practical on standard workstations. Each of the graphs comprising Fig. 3 required about two 
hours of CPU time on a large- memory (1Gb) DEC Alpha 500. The stability of Heisenberg 
spin ladder calculations in the absence of reorthogonalization contributes to the numerical 
efficiency: reorthogonalization can be costly because of the need for Lanczos vector storage 
and associated i/o operations. Thus such methods-of-moments calculations of dynamical 
spin responses for Heisenberg spin chains should be practical for L = 14 on workstations 
similar to ours, and extensions well beyond this possible with supercomputers. 

We thank J. Piekarewicz and J. Shepard for bringing this problem to our attention and 
for several helpful discussions. This work was supported in part by the US Department 
of Energy. One of us (DY) ackowledges the support of the University of Washington and 
National Science Foundation Physics Research Experiences for Undergraduates Program. 
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FIGURES 

FIG. 1. Results for 2 xL = 24 spin sites showing the convergence of the dynamical spin 
response per spin as a function of the number of iterations performed. The required number of 
iterations depends on the desired resolution, which in these calculations is determined by the choice 
of smearing function. A Gaussian with a = 0.05 has been used. 

FIG. 2. The dynamical spin response per spin for a ladder with 2 x L = 20 sites calculated 
with (70 iterations) and without (70 and 500 iterations) reorthogonalization. The results without 
reorthogonalization remain stable well past the point where S'(g7r,7r)'^) has fully converged. 

FIG. 3. The evolution of the spin response function per site for ladders with 2 x L = 24 sites 
as a function of J_\_/J\\- 



TABLES 

TABLE L The ground state energy per spin and interaction strength Eq/2LJ\\ for 2 x L 

Heisenberg spin ladders with J_|_ = Jii. 

L Present Barnes et al. 

I -0.6025112 -0.602511 

5 -0.5638793 

6 -0.5844372 -0.584437 

7 -0.5739430 

8 -0.5802030 -0.580203 

9 -0.5766331 

10 -0.5788595 -0.578860 

II -0.5775071 

12 -0.5783722 -0.578375 

13 -0.5778259 

14 -0.5781816 



TABLE II. As in Figure 1, only for the singlet-triplet gap {Ei — Eq)/J\\. 

L Present Barnes et al. 

I 0.8200894 0.820089 

5 0.8761249 

6 0.6265690 0.626570 

7 0.7734289 

8 0.5573976 0.557398 

9 0.7039126 

10 0.5281070 0.528106 

II 0.6558908 

12 0.5147836 0.514999 

13 0.6218955 

14 0.5084957 



TABLE III. The dynamical spin response per spin for the 2 x L = 24 Heisenberg spin ladder 
divided into the lowest triplet state contribution and the contribution carried by all higher states. 
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